source('../env.R')
community_data = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics.csv'))
Rows: 792 Columns: 93── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): is_urban_threshold
dbl (92): city_id, mntd_normalised, mntd_actual, mntd_min, mntd_max, mntd_mean, mntd_sd, fd_normalised, fd_actual, fd_min, fd_max, fd_mean, fd_sd, FRic_normal...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
community_data$is_urban_threshold = factor(community_data$is_urban_threshold, levels = c('low', 'medium', 'high'), labels = c('Low', 'Medium', 'High'))
head(community_data)
colnames(community_data)
[1] "city_id" "mntd_normalised" "mntd_actual" "mntd_min"
[5] "mntd_max" "mntd_mean" "mntd_sd" "fd_normalised"
[9] "fd_actual" "fd_min" "fd_max" "fd_mean"
[13] "fd_sd" "FRic_normalised" "FRic_actual" "FRic_min"
[17] "FRic_max" "FRic_mean" "FRic_sd" "mass_fd_normalised"
[21] "mass_fd_actual" "mass_fd_min" "mass_fd_max" "mass_fd_mean"
[25] "mass_fd_sd" "mass_var_normalised" "mass_var_actual" "mass_var_min"
[29] "mass_var_max" "mass_var_mean" "mass_var_sd" "mass_sdndr_normalised"
[33] "mass_sdndr_actual" "mass_sdndr_min" "mass_sdndr_max" "mass_sdndr_mean"
[37] "mass_sdndr_sd" "locomotory_trait_fd_normalised" "locomotory_trait_fd_actual" "locomotory_trait_fd_min"
[41] "locomotory_trait_fd_max" "locomotory_trait_fd_mean" "locomotory_trait_fd_sd" "locomotory_trait_var_normalised"
[45] "locomotory_trait_var_actual" "locomotory_trait_var_min" "locomotory_trait_var_max" "locomotory_trait_var_mean"
[49] "locomotory_trait_var_sd" "locomotory_trait_sdndr_normalised" "locomotory_trait_sdndr_actual" "locomotory_trait_sdndr_min"
[53] "locomotory_trait_sdndr_max" "locomotory_trait_sdndr_mean" "locomotory_trait_sdndr_sd" "trophic_trait_fd_normalised"
[57] "trophic_trait_fd_actual" "trophic_trait_fd_min" "trophic_trait_fd_max" "trophic_trait_fd_mean"
[61] "trophic_trait_fd_sd" "trophic_trait_var_normalised" "trophic_trait_var_actual" "trophic_trait_var_min"
[65] "trophic_trait_var_max" "trophic_trait_var_mean" "trophic_trait_var_sd" "trophic_trait_sdndr_normalised"
[69] "trophic_trait_sdndr_actual" "trophic_trait_sdndr_min" "trophic_trait_sdndr_max" "trophic_trait_sdndr_mean"
[73] "trophic_trait_sdndr_sd" "gape_width_fd_normalised" "gape_width_fd_actual" "gape_width_fd_min"
[77] "gape_width_fd_max" "gape_width_fd_mean" "gape_width_fd_sd" "gape_width_var_normalised"
[81] "gape_width_var_actual" "gape_width_var_min" "gape_width_var_max" "gape_width_var_mean"
[85] "gape_width_var_sd" "gape_width_sdndr_normalised" "gape_width_sdndr_actual" "gape_width_sdndr_min"
[89] "gape_width_sdndr_max" "gape_width_sdndr_mean" "gape_width_sdndr_sd" "urban_pool_size"
[93] "is_urban_threshold"
Join on realms
city_to_realm = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))
Rows: 342 Columns: 2── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
community_data_with_realm = left_join(community_data, city_to_realm)
Joining with `by = join_by(city_id)`
Cities as points
city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp')))
Warning: st_centroid assumes attributes are constant over geometries of x
city_points_low = city_points %>% left_join(community_data[community_data$is_urban_threshold == 'Low',])
Joining with `by = join_by(city_id)`
city_points_med = city_points %>% left_join(community_data[community_data$is_urban_threshold == 'Medium',])
Joining with `by = join_by(city_id)`
city_points_high = city_points %>% left_join(community_data[community_data$is_urban_threshold == 'High',])
Joining with `by = join_by(city_id)`
sf::sf_use_s2(FALSE)
Spherical geometry (s2) switched off
COUNTRY_BOUNDARIES = '/Users/james/Dropbox/PhD/WorldBank_countries_Admin0_10m/WB_countries_Admin0_10m.shp'
world_map = st_simplify(st_read(COUNTRY_BOUNDARIES), dTolerance = 0.02)
Reading layer `WB_countries_Admin0_10m' from data source `/Users/james/Dropbox/PhD/WorldBank_countries_Admin0_10m/WB_countries_Admin0_10m.shp' using driver `ESRI Shapefile'
Simple feature collection with 251 features and 52 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -180 ymin: -59.47275 xmax: 180 ymax: 83.6341
Geodetic CRS: WGS 84
Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees
normalised_colours_scale = scale_colour_gradient2(
low = "darkgreen",
mid = "yellow",
high = "red",
midpoint = 0.5,
space = "Lab",
na.value = "grey50",
guide = "colourbar",
aesthetics = "colour",
)
Load community data, and create long format version
communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))
Rows: 2462 Columns: 9── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, jetz_species_name, seasonal, presence, origin
dbl (1): city_id
lgl (3): present_urban_high, present_urban_med, present_urban_low
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
communities_summary = communities %>% group_by(city_id) %>% summarise(
regional_pool_size = n(),
urban_size_high = sum(present_urban_high),
urban_size_med = sum(present_urban_med),
urban_size_low = sum(present_urban_low)
)
communities_summary_long = bind_rows(
communities_summary %>% rename(urban_pool_size = 'urban_size_high') %>% dplyr::select(city_id, regional_pool_size, urban_pool_size) %>% mutate(is_urban_threshold = 'High'),
communities_summary %>% rename(urban_pool_size = 'urban_size_med') %>% dplyr::select(city_id, regional_pool_size, urban_pool_size) %>% mutate(is_urban_threshold = 'Medium'),
communities_summary %>% rename(urban_pool_size = 'urban_size_low') %>% dplyr::select(city_id, regional_pool_size, urban_pool_size) %>% mutate(is_urban_threshold = 'Low')
)
communities_summary_long
Load trait data
traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_jetz.csv'))
Rows: 304 Columns: 5── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): jetz_species_name
dbl (4): gape_width, trophic_trait, locomotory_trait, mass
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(traits)
Load realm geo
st_crs(resolve)
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
geom_normalised_histogram = function(name, gg, legend.position = "right") {
gg +
geom_histogram(aes(fill = core_realm), binwidth = 0.1, position = "dodge") +
geom_vline(aes(xintercept = 0.5), color = "#000000", size = 0.4) +
geom_vline(aes(xintercept = 0), color = "#000000", size = 0.2, linetype = "dashed") +
geom_vline(aes(xintercept = 1), color = "#000000", size = 0.2, linetype = "dashed") +
ylab("Number of cities") + xlab("Normalised Response") + ylim(c(0, 70)) +
labs(title = name, fill = 'Realm') +
facet_wrap(~ is_urban_threshold, ncol=1) +
theme_bw() +
theme(legend.position=legend.position)
}
geom_map = function(map_sf, title) {
norm_mntd_analysis_geo = ggplot() +
geom_sf(data = world_map, aes(geometry = geometry)) +
map_sf +
normalised_colours_scale +
labs(title = title, colour = 'Normalised\nResponse')
}
norm_mntd_analysis_plot = geom_normalised_histogram(
'MNTD',
ggplot(community_data_with_realm, aes(mntd_normalised))
)
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.
norm_mntd_analysis_plot
norm_mntd_analysis_geo = geom_map(geom_sf(data = city_points_low, aes(color = mntd_normalised, geometry = geometry)), 'Low threshold')
norm_mntd_analysis_geo
ggarrange(norm_mntd_analysis_plot,
ggarrange(
geom_normalised_histogram('MNTD', ggplot(community_data_with_realm, aes(mntd_normalised)), legend.position = "bottom"),
geom_map(geom_sf(data = city_points_med, aes(color = mntd_normalised, geometry = geometry)), 'Medium threshold'),
geom_map(geom_sf(data = city_points_high, aes(color = mntd_normalised, geometry = geometry)), 'High threshold'),
ncol=1, common.legend = T),
ncol = 2)
ggsave(filename(FIGURES_OUTPUT_DIR, 'normalised_mntd.jpg'), width = 2500, height=2500, units = 'px')
norm_fd_analysis_plot = geom_normalised_histogram(
'FD',
ggplot(community_data_with_realm, aes(fd_normalised))
)
norm_fd_analysis_plot
norm_fd_analysis_geo = geom_map(geom_sf(data = city_points_low, aes(color = fd_normalised, geometry = geometry)), 'Low threshold')
norm_fd_analysis_geo
ggarrange(geom_normalised_histogram('FD', ggplot(community_data_with_realm, aes(fd_normalised)), legend.position = "bottom"),
ggarrange(
norm_fd_analysis_geo,
geom_map(geom_sf(data = city_points_med, aes(color = fd_normalised, geometry = geometry)), 'Medium threshold'),
geom_map(geom_sf(data = city_points_high, aes(color = fd_normalised, geometry = geometry)), 'High threshold'),
ncol=1, common.legend = T),
ncol = 2)
ggsave(filename(FIGURES_OUTPUT_DIR, 'normalised_fd.jpg'), width = 2500, height=2500, units = 'px')
norm_fd_loco_analysis_plot = geom_normalised_histogram(
'FD Locomotory',
ggplot(community_data_with_realm, aes(locomotory_trait_fd_normalised))
)
norm_fd_loco_analysis_plot
norm_var_loco_analysis_plot = geom_normalised_histogram(
'Var Locomotory',
ggplot(community_data_with_realm, aes(locomotory_trait_var_normalised))
)
norm_var_loco_analysis_plot
norm_fd_trophic_analysis_plot = geom_normalised_histogram(
'FD Trophic',
ggplot(community_data_with_realm, aes(trophic_trait_fd_normalised))
)
norm_fd_trophic_analysis_plot
norm_var_trophic_analysis_plot = geom_normalised_histogram(
'Var Trophic',
ggplot(community_data_with_realm, aes(trophic_trait_var_normalised))
)
norm_var_trophic_analysis_plot
norm_fd_gape_analysis_plot = geom_normalised_histogram(
'FD Gape Width',
ggplot(community_data_with_realm, aes(gape_width_fd_normalised))
)
norm_fd_gape_analysis_plot
norm_var_gape_analysis_plot = geom_normalised_histogram(
'Var Gape Width',
ggplot(community_data_with_realm, aes(gape_width_var_normalised))
)
norm_var_gape_analysis_plot
norm_fd_mass_analysis_plot = geom_normalised_histogram(
'FD Mass',
ggplot(community_data_with_realm, aes(mass_fd_normalised))
)
norm_fd_mass_analysis_plot
norm_var_mass_analysis_plot = geom_normalised_histogram(
'Var Mass',
ggplot(community_data_with_realm, aes(mass_var_normalised))
)
norm_var_mass_analysis_plot
ggplot(community_data_with_realm, aes(x = fd_normalised, y = mntd_normalised, colour = core_realm)) +
geom_point() +
ylab("MNTD") +
xlab("FD") +
theme_bw() +
facet_wrap(~ is_urban_threshold, ncol = 2) + labs(color = "Realm")
ggsave(filename(FIGURES_OUTPUT_DIR, 'mntd_by_fd.jpg'))
Saving 7.29 x 4.51 in image
mntd_fd_analysis = community_data_with_realm %>%
dplyr::select(city_id, is_urban_threshold, mntd_normalised, fd_normalised) %>%
left_join(communities_summary_long) %>%
mutate(urban_pool_perc = urban_pool_size * 100 / regional_pool_size)
Joining with `by = join_by(city_id, is_urban_threshold)`
mntd_fd_analysis
ggpairs(mntd_fd_analysis %>% dplyr::filter(is_urban_threshold == 'Low') %>% dplyr::select(mntd_normalised, fd_normalised, regional_pool_size, urban_pool_size, urban_pool_perc))
ggpairs(mntd_fd_analysis %>% dplyr::filter(is_urban_threshold == 'Medium') %>% dplyr::select(mntd_normalised, fd_normalised, regional_pool_size, urban_pool_size, urban_pool_perc))
ggsave(filename(FIGURES_OUTPUT_DIR, 'mntd_by_fd_pairs_medium.jpg'))
Saving 7.29 x 4.51 in image
ggpairs(mntd_fd_analysis %>% dplyr::filter(is_urban_threshold == 'High') %>% dplyr::select(mntd_normalised, fd_normalised, regional_pool_size, urban_pool_size, urban_pool_perc))
communities_with_traits = communities %>% left_join(traits) %>% left_join(city_to_realm) %>% filter(core_realm != 'Oceania')
Joining with `by = join_by(jetz_species_name)`Joining with `by = join_by(city_id)`
head(communities_with_traits)
realms = unique(communities_with_traits$core_realm)
realms
[1] "Nearctic" "Neotropic" "Palearctic" "Afrotropic" "Indomalayan" "Australasia"
communities_with_traits_long = bind_rows(
communities_with_traits %>% dplyr::select(gape_width, trophic_trait, locomotory_trait, mass, core_realm) %>% mutate(record_type = 'Regional'),
communities_with_traits %>% filter(present_urban_low) %>% dplyr::select(gape_width, trophic_trait, locomotory_trait, mass, core_realm) %>% mutate(record_type = 'Urban Low Threshold'),
communities_with_traits %>% filter(present_urban_med) %>% dplyr::select(gape_width, trophic_trait, locomotory_trait, mass, core_realm) %>% mutate(record_type = 'Urban Medium Threshold'),
communities_with_traits %>% filter(present_urban_high) %>% dplyr::select(gape_width, trophic_trait, locomotory_trait, mass, core_realm) %>% mutate(record_type = 'Urban High Threshold')
)
communities_with_traits_long$record_type = factor(communities_with_traits_long$record_type, levels = c('Regional', 'Urban Low Threshold', 'Urban Medium Threshold', 'Urban High Threshold'))
head(communities_with_traits_long)
convex_hull_loco_trophic_per_realm = function(filtered_df, realms) {
result = data.frame()
for (realm in realms) {
result = rbind(result,
filtered_df %>%
filter(core_realm == realm) %>%
slice(chull(trophic_trait, locomotory_trait)) %>%
dplyr::select(trophic_trait, locomotory_trait) %>%
mutate(core_realm = realm)
)
}
result
}
regional_hull = convex_hull_loco_trophic_per_realm(communities_with_traits, realms)
urban_hull_high = convex_hull_loco_trophic_per_realm(communities_with_traits %>% filter(present_urban_high), realms)
urban_hull_med = convex_hull_loco_trophic_per_realm(communities_with_traits %>% filter(present_urban_med), realms)
urban_hull_low = convex_hull_loco_trophic_per_realm(communities_with_traits %>% filter(present_urban_low), realms)
ggplot(data = communities_with_traits_long, aes(x = trophic_trait, y = locomotory_trait)) +
geom_polygon(data = regional_hull, alpha = 0.1, fill = "green", color="green") +
geom_polygon(data = urban_hull_low, alpha = 0.15, fill = "yellow", color="yellow") +
geom_polygon(data = urban_hull_med, alpha = 0.2, fill = "orange", color = "orange") +
geom_polygon(data = urban_hull_high, alpha = 0.25, fill = "red", color = "red") +
geom_point(aes(colour = record_type)) +
theme_bw() + scale_color_manual(values=c("green", "yellow", "orange", "red")) +
xlab('Trophic Trait') + ylab('Locomotory Trait') +
theme(legend.position="bottom") +
facet_wrap(~ core_realm)
ggsave(filename(FIGURES_OUTPUT_DIR, 'traits_by_realm_loco_trophic.jpg'))
Saving 7.29 x 4.51 in image
convex_hull_gape_mass_per_realm = function(filtered_df, realms) {
result = data.frame()
for (realm in realms) {
result = rbind(result,
filtered_df %>%
filter(core_realm == realm) %>%
slice(chull(gape_width, mass)) %>%
dplyr::select(gape_width, mass) %>%
mutate(core_realm = realm)
)
}
result
}
regional_hull_gm = convex_hull_gape_mass_per_realm(communities_with_traits, realms)
urban_hull_high_gm = convex_hull_gape_mass_per_realm(communities_with_traits %>% filter(present_urban_high), realms)
urban_hull_med_gm = convex_hull_gape_mass_per_realm(communities_with_traits %>% filter(present_urban_med), realms)
urban_hull_low_gm = convex_hull_gape_mass_per_realm(communities_with_traits %>% filter(present_urban_low), realms)
ggplot(data = communities_with_traits_long, aes(x = gape_width, y = mass)) +
geom_polygon(data = regional_hull_gm, alpha = 0.1, fill = "green", color="green") +
geom_polygon(data = urban_hull_low_gm, alpha = 0.15, fill = "yellow", color="yellow") +
geom_polygon(data = urban_hull_med_gm, alpha = 0.2, fill = "orange", color = "orange") +
geom_polygon(data = urban_hull_high_gm, alpha = 0.25, fill = "red", color = "red") +
geom_point(aes(colour = record_type)) +
theme_bw() + scale_color_manual(values=c("green", "yellow", "orange", "red")) +
xlab('Gape Width') + ylab('Mass') +
theme(legend.position="bottom") +
facet_wrap(~ core_realm)
ggsave(filename(FIGURES_OUTPUT_DIR, 'traits_by_realm_gape_mass.jpg'))
Saving 7.29 x 4.51 in image
phylo_tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
ggtree(phylo_tree, layout='circular')
plot_phylo = function(realm, community_df = communities_with_traits) {
species_df = community_df %>% filter(core_realm == realm) %>% group_by(jetz_species_name) %>% summarise(
regional_pools = n(),
urban_pools_high = sum(present_urban_high),
urban_pools_medium = sum(present_urban_med),
urban_pools_low = sum(present_urban_low)
) %>% mutate(
urban_tolerence_high = urban_pools_high / regional_pools,
urban_tolerence_med = urban_pools_medium / regional_pools,
urban_tolerence_low = urban_pools_low / regional_pools
)
tree_cropped <- ladderize(drop.tip(phylo_tree, setdiff(phylo_tree$tip.label, species_df$jetz_species_name)))
p = ggtree(tree_cropped) + geom_tiplab(align=T)
p1 = facet_plot(p, panel='High', data=species_df, geom=geom_segment, aes(x=0, xend=urban_tolerence_high, y=y, yend=y), size=3, color='red')
p2 = facet_plot(p1, panel='Medium', data=species_df, geom=geom_segment, aes(x=0, xend=urban_tolerence_med, y=y, yend=y), size=3, color='orange')
p3 = facet_plot(p2, panel='Low', data=species_df, geom=geom_segment, aes(x=0, xend=urban_tolerence_low, y=y, yend=y), size=3, color='yellow')
facet_widths(p3 + xlim_tree(60) + theme_tree2(), c(Tree = 5)) + labs(title = realm, subtitle = 'Urban tolerance')
}
realms
[1] "Nearctic" "Neotropic" "Palearctic" "Afrotropic" "Indomalayan" "Australasia"
plot_phylo('Nearctic')
ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny_nearctic.jpg'))
Saving 7.29 x 4.51 in image
plot_phylo('Neotropic')
ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny_neotropic.jpg'), width=2187, units="px")
Saving 2187 x 3000 px image
plot_phylo('Palearctic')
ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny_palearctic.jpg'))
Saving 7.29 x 4.51 in image
plot_phylo('Afrotropic')
ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny_afrotropic.jpg'))
Saving 7.29 x 4.51 in image
plot_phylo('Australasia')
ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny_australasia.jpg'))
Saving 7.29 x 4.51 in image
It seems reasonable to expect that cities with simialr regional pools will have similar species entering the city, and thus a similar response to urbanisation.
to_species_matrix = function(filtered_communities) {
filtered_communities %>%
dplyr::select(city_id, jetz_species_name) %>%
distinct() %>%
mutate(present = TRUE) %>%
pivot_wider(
names_from = jetz_species_name,
values_from = "present",
values_fill = list(present = F)
) %>%
tibble::column_to_rownames(var='city_id')
}
community_nmds = function(filtered_communities) {
species_matrix = to_species_matrix(filtered_communities)
nmds <- metaMDS(species_matrix, k=2, trymax = 30)
nmds_result = data.frame(scores(nmds))
nmds_result$city_id = as.double(rownames(scores(nmds)))
rownames(nmds_result) = NULL
nmds_result %>%
left_join(city_points[,c('city_id', 'city_nm')])
}
plot_nmds = function(plot_data, clusters, realm, plot_is_urban_threshold) {
pd = plot_data %>%
left_join(community_data) %>%
filter(is_urban_threshold == plot_is_urban_threshold)
ggplot(pd, aes(x = NMDS1, y = NMDS2)) +
geom_polygon(data = clusters, aes(fill = cluster), alpha = 0.2) +
geom_jitter(aes(colour = mntd_normalised, size = fd_normalised)) +
normalised_colours_scale +
geom_text_repel(aes(label = city_nm), max.overlaps = 25) +
labs(
title = realm,
subtitle = paste('Regional pool, with urban community values using', plot_is_urban_threshold, 'threshold'),
color = 'MNTD Normalised',
size = 'FD Normalised')
}
https://www.datacamp.com/tutorial/k-means-clustering-r
scree_plot = function(community_nmds_data) {
# Decide how many clusters to look at
n_clusters <- 10
# Initialize total within sum of squares error: wss
wss <- numeric(n_clusters)
set.seed(123)
# Look over 1 to n possible clusters
for (i in 1:n_clusters) {
# Fit the model: km.out
km.out <- kmeans(community_nmds_data[,c('NMDS1','NMDS2')], centers = i, nstart = 20)
# Save the within cluster sum of squares
wss[i] <- km.out$tot.withinss
}
# Produce a scree plot
wss_df <- tibble(clusters = 1:n_clusters, wss = wss)
scree_plot <- ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
geom_point(size = 4)+
geom_line() +
scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
xlab('Number of clusters')
scree_plot
}
convex_hull_by_cluster = function(community_df) {
result = data.frame()
clusters = unique(community_df$cluster)
for (cluster_i in clusters) {
result = rbind(result,
community_df %>%
filter(cluster == cluster_i) %>%
slice(chull(NMDS1, NMDS2)) %>%
dplyr::select(NMDS1, NMDS2, cluster)
)
}
result
}
nearctic_communities = community_nmds(communities_with_traits %>% filter(core_realm == 'Nearctic'))
Run 0 stress 0.06772313
Run 1 stress 0.09925753
Run 2 stress 0.0814111
Run 3 stress 0.0714024
Run 4 stress 0.08477133
Run 5 stress 0.06792384
... Procrustes: rmse 0.002677399 max resid 0.01592546
Run 6 stress 0.07140237
Run 7 stress 0.09410676
Run 8 stress 0.06772314
... Procrustes: rmse 0.000007756193 max resid 0.00001712457
... Similar to previous best
Run 9 stress 0.0847715
Run 10 stress 0.06772321
... Procrustes: rmse 0.00005153412 max resid 0.00008951089
... Similar to previous best
Run 11 stress 0.06792388
... Procrustes: rmse 0.00267859 max resid 0.01598036
Run 12 stress 0.0677232
... Procrustes: rmse 0.00004962282 max resid 0.00008595323
... Similar to previous best
Run 13 stress 0.08193052
Run 14 stress 0.09279085
Run 15 stress 0.08494123
Run 16 stress 0.09422939
Run 17 stress 0.07122361
Run 18 stress 0.08521544
Run 19 stress 0.07427005
Run 20 stress 0.06792383
... Procrustes: rmse 0.002676712 max resid 0.01589483
*** Solution reached
Joining with `by = join_by(city_id)`
nearctic_communities
scree_plot(nearctic_communities)
Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3300)Warning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3300)
nearctic_clusters <- kmeans(nearctic_communities[,c('NMDS1', 'NMDS2')], centers = 3, nstart = 20)
nearctic_communities$cluster = as.factor(nearctic_clusters$cluster)
ggplot(nearctic_communities, aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point()
nearctic_cluster_hulls = convex_hull_by_cluster(nearctic_communities)
plot_nmds(nearctic_communities, nearctic_cluster_hulls, 'Nearctic', 'High')
Joining with `by = join_by(city_id)`
plot_nmds(nearctic_communities, nearctic_cluster_hulls, 'Nearctic', 'Medium')
Joining with `by = join_by(city_id)`
plot_nmds(nearctic_communities, nearctic_cluster_hulls, 'Nearctic', 'Low')
Joining with `by = join_by(city_id)`
species_in_neartic = nearctic_communities %>% dplyr::select(city_id, cluster) %>% left_join(communities_with_traits)
Joining with `by = join_by(city_id)`
plot_phylo('Nearctic', species_in_neartic %>% filter(cluster == 1))
ggplot() +
geom_sf(data = resolve_nearctic, aes(geometry = geometry), alpha = 0.2) +
geom_sf(data = nearctic_communities %>% filter(cluster == 1) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
normalised_colours_scale + theme_bw()
Joining with `by = join_by(city_id, city_nm, geometry)`
plot_phylo('Nearctic', species_in_neartic %>% filter(cluster == 2))
ggplot() +
geom_sf(data = resolve_nearctic, aes(geometry = geometry), alpha = 0.2) +
geom_sf(data = nearctic_communities %>% filter(cluster == 2) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
normalised_colours_scale + theme_bw()
Joining with `by = join_by(city_id, city_nm, geometry)`
ggplot() +
geom_sf(data = resolve_nearctic, aes(geometry = geometry), alpha = 0.2) +
geom_sf(data = nearctic_communities %>% filter(cluster == 3) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
normalised_colours_scale + theme_bw()
Joining with `by = join_by(city_id, city_nm, geometry)`
neotropic_communities = community_nmds(communities_with_traits %>% filter(core_realm == 'Neotropic'))
Run 0 stress 0.134619
Run 1 stress 0.1344509
... New best solution
... Procrustes: rmse 0.006928659 max resid 0.0457892
Run 2 stress 0.1369387
Run 3 stress 0.134619
... Procrustes: rmse 0.006940788 max resid 0.04575861
Run 4 stress 0.1344509
... New best solution
... Procrustes: rmse 0.00004403156 max resid 0.0001387448
... Similar to previous best
Run 5 stress 0.1346406
... Procrustes: rmse 0.007192008 max resid 0.0457332
Run 6 stress 0.1346226
... Procrustes: rmse 0.007280471 max resid 0.04573126
Run 7 stress 0.136377
Run 8 stress 0.1348046
... Procrustes: rmse 0.006663518 max resid 0.04610571
Run 9 stress 0.1406945
Run 10 stress 0.1346226
... Procrustes: rmse 0.007291867 max resid 0.04571052
Run 11 stress 0.134619
... Procrustes: rmse 0.006932652 max resid 0.04571971
Run 12 stress 0.134433
... New best solution
... Procrustes: rmse 0.001186888 max resid 0.006991429
... Similar to previous best
Run 13 stress 0.1348237
... Procrustes: rmse 0.006651112 max resid 0.04607972
Run 14 stress 0.1346405
... Procrustes: rmse 0.007243159 max resid 0.04569399
Run 15 stress 0.1406945
Run 16 stress 0.134619
... Procrustes: rmse 0.006830928 max resid 0.04574414
Run 17 stress 0.1344509
... Procrustes: rmse 0.001186553 max resid 0.00692233
... Similar to previous best
Run 18 stress 0.141222
Run 19 stress 0.1346207
... Procrustes: rmse 0.006999734 max resid 0.04555
Run 20 stress 0.134433
... New best solution
... Procrustes: rmse 0.00002505145 max resid 0.00008159977
... Similar to previous best
*** Solution reached
Joining with `by = join_by(city_id)`
neotropic_communities
scree_plot(neotropic_communities)
neotropic_clusters <- kmeans(neotropic_communities[,c('NMDS1', 'NMDS2')], centers = 3, nstart = 20)
neotropic_communities$cluster = as.factor(neotropic_clusters$cluster)
ggplot(neotropic_communities, aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point()
neotropic_cluster_hulls = convex_hull_by_cluster(neotropic_communities)
plot_nmds(neotropic_communities, neotropic_cluster_hulls, 'Neotropic', 'High')
Joining with `by = join_by(city_id)`
plot_nmds(neotropic_communities, neotropic_cluster_hulls, 'Neotropic', 'Medium')
Joining with `by = join_by(city_id)`
plot_nmds(neotropic_communities, neotropic_cluster_hulls, 'Neotropic', 'Low')
Joining with `by = join_by(city_id)`
species_in_neotropic = neotropic_communities %>% dplyr::select(city_id, cluster) %>% left_join(communities_with_traits)
Joining with `by = join_by(city_id)`
plot_phylo('Neotropic', species_in_neotropic %>% filter(cluster == 1))
resolve_neotropic = resolve %>% filter(REALM == 'Neotropic')
ggplot() +
geom_sf(data = resolve_neotropic, aes(geometry = geometry), alpha = 0.2) +
geom_sf(data = neotropic_communities %>% filter(cluster == 1) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
normalised_colours_scale + theme_bw()
Joining with `by = join_by(city_id, city_nm, geometry)`
plot_phylo('Neotropic', species_in_neotropic %>% filter(cluster == 2))
ggplot() +
geom_sf(data = resolve_neotropic, aes(geometry = geometry), alpha = 0.2) +
geom_sf(data = neotropic_communities %>% filter(cluster == 2) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
normalised_colours_scale + theme_bw()
Joining with `by = join_by(city_id, city_nm, geometry)`
plot_phylo('Neotropic', species_in_neotropic %>% filter(cluster == 3))
ggplot() +
geom_sf(data = resolve_neotropic, aes(geometry = geometry), alpha = 0.2) +
geom_sf(data = neotropic_communities %>% filter(cluster == 3) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
normalised_colours_scale + theme_bw()
Joining with `by = join_by(city_id, city_nm, geometry)`
indomalayan_communities = community_nmds(communities_with_traits %>% filter(core_realm == 'Indomalayan'))
Run 0 stress 0.1193009
Run 1 stress 0.1285877
Run 2 stress 0.1472786
Run 3 stress 0.1289009
Run 4 stress 0.1240017
Run 5 stress 0.157254
Run 6 stress 0.1509529
Run 7 stress 0.161852
Run 8 stress 0.1558148
Run 9 stress 0.1239373
Run 10 stress 0.1165214
... New best solution
... Procrustes: rmse 0.03581633 max resid 0.2499892
Run 11 stress 0.1151087
... New best solution
... Procrustes: rmse 0.02627242 max resid 0.2457618
Run 12 stress 0.1167036
Run 13 stress 0.1220414
Run 14 stress 0.1229737
Run 15 stress 0.1382712
Run 16 stress 0.1472537
Run 17 stress 0.1502501
Run 18 stress 0.1666794
Run 19 stress 0.1398437
Run 20 stress 0.1448979
Run 21 stress 0.1252406
Run 22 stress 0.1443181
Run 23 stress 0.1156951
Run 24 stress 0.1618047
Run 25 stress 0.1512284
Run 26 stress 0.1245209
Run 27 stress 0.1166941
Run 28 stress 0.1156757
Run 29 stress 0.1257237
Run 30 stress 0.1169841
*** No convergence -- monoMDS stopping criteria:
1: no. of iterations >= maxit
26: stress ratio > sratmax
3: scale factor of the gradient < sfgrmin
Joining with `by = join_by(city_id)`
indomalayan_communities
scree_plot(indomalayan_communities)
indomalayan_clusters <- kmeans(indomalayan_communities[,c('NMDS1', 'NMDS2')], centers = 4, nstart = 20)
indomalayan_communities$cluster = as.factor(indomalayan_clusters$cluster)
ggplot(indomalayan_communities, aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point()
indomalayan_cluster_hulls = convex_hull_by_cluster(indomalayan_communities)
plot_nmds(indomalayan_communities, indomalayan_cluster_hulls, 'Indomalayan', 'High')
Joining with `by = join_by(city_id)`
plot_nmds(indomalayan_communities, indomalayan_cluster_hulls, 'Indomalayan', 'Medium')
Joining with `by = join_by(city_id)`
plot_nmds(indomalayan_communities, indomalayan_cluster_hulls, 'Indomalayan', 'Low')
Joining with `by = join_by(city_id)`
species_in_indomalayan = indomalayan_communities %>% dplyr::select(city_id, cluster) %>% left_join(communities_with_traits)
Joining with `by = join_by(city_id)`
plot_phylo('Indomalayan', species_in_indomalayan %>% filter(cluster == 1))
resolve_indomalayan = resolve %>% filter(REALM == 'Indomalayan')
ggplot() +
geom_sf(data = resolve_indomalayan, aes(geometry = geometry), alpha = 0.2) +
geom_sf(data = indomalayan_communities %>% filter(cluster == 1) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
normalised_colours_scale + theme_bw()
Joining with `by = join_by(city_id, city_nm, geometry)`
ggplot() +
geom_sf(data = resolve_indomalayan, aes(geometry = geometry), alpha = 0.2) +
geom_sf(data = indomalayan_communities %>% filter(cluster == 2) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
normalised_colours_scale + theme_bw()
Joining with `by = join_by(city_id, city_nm, geometry)`
plot_phylo('Indomalayan', species_in_indomalayan %>% filter(cluster == 3))
ggplot() +
geom_sf(data = resolve_indomalayan, aes(geometry = geometry), alpha = 0.2) +
geom_sf(data = indomalayan_communities %>% filter(cluster == 3) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
normalised_colours_scale + theme_bw()
Joining with `by = join_by(city_id, city_nm, geometry)`
plot_phylo('Indomalayan', species_in_indomalayan %>% filter(cluster == 4))
ggplot() +
geom_sf(data = resolve_indomalayan, aes(geometry = geometry), alpha = 0.2) +
geom_sf(data = indomalayan_communities %>% filter(cluster == 4) %>% left_join(city_points_med), aes(color = mntd_normalised, geometry = geometry)) +
normalised_colours_scale + theme_bw()
Joining with `by = join_by(city_id, city_nm, geometry)`